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1.  Introduction 


Future  military  analog  fiber-optic  links  require  high-bandwidth  (20-30  GHz)  and  large 
dynamic  range  (>120  dB)  that  are,  furthermore,  compact,  light  weight,  and  with  low 
power  consumption.  Technologies  that  have  been  investigated  so  far  can  not  meet  the 
above-mentioned  demands.  These  technologies  include,  cw  laser  diode  direct 
modulation,  mode-locked  laser  diode  external  modulation,  and  external  modulation. 

Optical  coupling  of  a  Laser  diode  (LD)  with  an  external  cavity  (EC),  such  as 
shown  schematically  in  Fig.  1,  has  been  used  to  improve  laser  frequency  stability  and  to 
narrow  laser  linewidth.  Such  coupling  can  also  significantly  change  LD  dynamics.  For 
long  (>1  cm)  external  cavity  LDs  these  changes  lead  to  chaotic  self-modulation  of  laser 
light  with  dramatic  broadening  of  laser  spectrum  ("coherence-collapse"  effect).  For  a 
short  EC  LD  there  exists  another  type  of  instability  associated  with  the  EC  mode 
interactions  alone  [1].  The  compound-cavity  mode  interaction  inside  the  regime  of 
stable  CW  operation  lead  to  the  generation  of  high  frequency  (>20  GHz)  microwave 
oscillations.  Such  high-frequency  oscillations  arise,  as  has  been  conjectured  [1],  if  the 
cavity  losses  for  the  main  and  side  compound-cavity  modes  of  the  short  EC  LD  become 
equal  with  increasing  feedback.  This  equalization  does  not  lead,  however,  to  a  switching 
of  lasing  to  this  sidemode,  but  instead  induces  a  self-mode-locking  regime  with  deep 
oscillations  of  the  laser  output. 
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We  propose  to  use  short  EC  LDs  as  tunable  microwave  oscillators  with  frequency 
modulation  capabilities.  There  are  two  ways  the  frequency  can  be  modulated  [2]. 

•  It  has  been  shown  that  the  oscillation  frequency  increases  with  increasing 
feedback.  Therefore,  the  frequency  of  the  microwave  oscillation  can  be 
modulated  by  the  feedback 

•  It  should  be  noted  that  the  oscillation  is  a  result  of  the  interaction  of  compound- 
cavity  mode  interaction.  Therefore,  once  the  optical  length  of  the  individual 
cavity  is  changed,  the  oscillation  is  tuned. 

We  have  previously  computed  the  dynamics  of  short  EC  LDs,  mapped  out  the  system 
parameter  space  for  microwave  oscillations,  and  determined  the  frequency  range  of  the 
microwave  oscillations. 

The  ultimate  research  objective  is  to  experimentally  demonstrate  the  generation  of 
tunable  microwave  intensity  oscillations  in  laser  diodes  with  short  external  cavity  for  use 
in  analog  fiber-optic  links.  A  better  understanding  of  the  fundamental  mechanism  is 
therefore  essential  for  progress  towards  employing  these  devices  in  real  applications.  In 
this  report,  we  document  enhanced  understanding  of  the  underlying  physical  principle 
for  the  microwave-frequency  intensity  modulation  from  a  nonlinear  dynamics  point  of 
view,  and  ways  to  experimentally  measure  the  noise  characteristics  of  the  oscillations. 
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2.  Investigation  of  the  Noise  Characteristics 

2.1  Preliminary  Comment 

Periodic  orbits  in  nonlinear  optical  systems,  such  as  the  high-frequency  oscillations 
described  in  Sections  2  and  3  can  serve  as  subcarriers  in  optical  communications. 
Subcarrier  optical  communication  using  one  form  of  periodic  orbits  commonly  known  as 
self- sustained  pulsations  (SSP)  [5]  has  been  experimentally  demonstrated  both  for 
analog  [6]  and  digital  [7]  applications.  The  signal-to-noise  ratio  (bit-error-rate)  of  the 
analog  (digital)  subcarrier  systems  depends  on  the  nature  and  magnitude  of  the  noise  of 
the  subcarrier.  Since  only  one  harmonic  in  SSP  is  used  as  the  subcarrier,  the  amplitude 
and  phase  noise  of  each  harmonic  in  SSP  needs  to  be  characterized.  In  this  report  we 
report  the  first,  to  the  best  of  our  knowledge,  calibrated  measurement  of  amplitude  and 
phase  noise  of  self-pulsations  in  two-section  Fabry-Perot  laser  diodes  (LDs). 

2.2  Experimental  Setup 

Figure  1  illustrates  the  experimental  setup  for  the  measurement  of  amplitude  and  phase 
noise  of  SSP  in  commercially  available  780  nm  SSP  LDs.  The  SSP  is  first  converted  into 
RF  signals  with  a  high-speed  photodetector  (3  dB  bandwidth  >14  GHz).  The  RF 
spectrum  of  SSP  is  monitored  using  a  broadband  RF  spectrum  analyzer  (HP  8565E).  By 
using  bandpass  filters,  individual  harmonic  of  the  RF  signal  can  be  selected  for 
amplitude  and  phase  noise  measurement  using  the  HP  3048A  Noise  Measurement 
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System,  which  is  a  calibrated  measurement  system  dedicated  for  RF/microwave  oscillator 
noise  characterization  for  frequency  offsets  as  small  as  1  Hz. 


Figure.  1  Experimental  setup  for  the  amplitude  and  phase  noise  measurement 


2.3  Results 


RF  Spectrum 

986  990  994  998  1002  1006 


Frequency  (GHz) 


Figure.  2  RF  spectrum  of  SSP  with  a  1  GHz  fundamental  frequency  in  a  laser  diode. 
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A  typical  RF  spectrum  of  SSP  is  shown  in  figure  2.  At  70  mA  DC  bias  current,  the  SSP 
frequency  is  around  1  GHz  (The  threshold  is  at  50  mA).  More  than  15  harmonics  can 
be  observed.  Each  harmonic  has  a  certain  linewidth  which  increases  with  its  order. 
Details  around  the  fundamental  harmonic,  given  in  the  inset,  reveals  a  3-dB  linewidth  on 
the  order  of  1  MHz.  Because  the  fundamental  harmonic  was  used  in  previous 
experiments  [6,7]  and  because  of  the  instrumentation  frequency  limitations,  we  focused 
on  the  noise  characterization  of  the  fundamental  harmonic  in  the  present  effort.  With  a 
1  GHz  bandpass  filter  of  400  MHz  bandwidth,  the  amplitude  and  phase  noise  of 
fundamental  harmonic  are  investigated  nearly  independent  of  each. 


Frequency  Offset  (Hz) 

Figure.  3  Amplitude  noise  spectrum  of  the  fundamental  harmonic. 


The  amplitude  noise  is  measured  first  because  the  AM  noise  interface  offers  a 
better  than  110-dB  discrimination  against  phase  noise.  The  amplitude  noise  spectral 
density  in  dBc,  which  represents  the  AM  noise  sideband  to  carrier  power  ratio,  is  given  in 
figure  3.  Comparing  the  RF  and  AM  noise  spectra,  it  is  apparent  that  amplitude  noise 
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has  a  negligible  contribution  to  the  overall  noise  power  leading  to  the  1  MHz  linewidth. 
Therefore,  the  noise  of  the  fundamental  harmonic  has  to  be  primarily  phase  noise.  The 
phase  noise  interface  is  essentially  a  delay-line  frequency  discriminator,  which  measures 
the  noise  frequency  deviation  from  the  center  frequency  due  to  phase  noise.  Figure  4 
shows  the  phase  noise  spectrum  of  the  fundamental  harmonic  for  1  Hz  to  40  MHz 
frequency  offset  from  the  center  frequency  of  the  carrier.  Below  1  kHz  offset,  the  noise 
frequency  deviation  is  1/f  in  nature.  The  peak  around  60  Hz  results  from  the  power 
supply.  For  the  frequency  offsets  above  10  kHz,  frequency  deviation  noise  is  white  with 
a  magnitude  of  around  600  Hz  /  ~J Hz  .  The  decrease  near  40  MHz  offsets  is  due  to  the 
bandwidth  limitation  of  the  HP  3048A  system.  Assuming  a  600  Hz  1 4Hz  pure 
frequency  fluctuation  which  is  white  in  nature,  the  theoretically  [8]  estimated  carrier 
linewidth  is  about  1.13  MHz  which  matches  very  well  with  the  observed  KF  spectrum. 


Figure.  4  Phase  noise  spectrum  of  the  fundamental  harmonic 
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2.4  Summary 


We  have  investigated  experimentally  the  noise  behavior  of  self- sustained  pulsation  in 
laser  diodes.  Frequency  deviation  and  amplitude  fluctuation  of  the  pulsation  have  been 
characterized  in  the  range  of  1  Hz  to  40  MHz.  Even  though  phase  noise  dominates  over 
amplitude  noise,  both  are  1/f  in  nature  for  offsets  below  1  kHz,  and  white  for  offsets 
above  10  kHz.  This  technique  can  be  applied  to  the  noise  measurement  of  the  high- 
frequency  oscillation  in  ECLDs  also..  In  addition  to  understanding  and  ultimately 
improving  the  performance  of  periodic-orbit  subcarrier  optical  communication  systems, 
further  research  along  this  line  will  be  useful  for  other  applications  such  as  all-optical 
clock  recovery  [9]. 


3.  Investigation  of  Deterministic  Properties:  Formulation 
3.1  Rate  Equation 


Laser  Diode  cavity  External  Cavity 


R1  R 2  R3 


L  Lext 


Fig.  5  Schematic  of  a  laser  diode  with  an  external  cavity. 
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The  schematic  of  a  semiconductor  laser  diode  with  an  external  cavity  is  shown  in 
Fig.  1.  Rl  and  R2  are  the  reflectivity  of  the  laser  facets  and  R3  is  the  reflectivity  of  the 
external  mirror.  L  and  Lext  are  the  lengths  of  the  laser  diode  cavity  and  the  EC;  fext  is 
the  feedback.  We  will  consider  <<1,  where  t  —  2  Lex{  I  c  is  the  photon  round-trip 
time  in  EC  and  co  R  is  the  angular  frequency  of  the  relaxation  oscillation. 


We  adopt  the  Lang  and  Kobayashi  rate  equations  for  the  SEC  LD.  There  are 
several  assumptions  for  these  rate  equations.  First,  the  solitary  laser  is  assumed  to  be 
single  mode,  which  implies  r»  r0  (  r0  =  2  Lng  /  c  ,  is  the  photon  round-trip  time  in  the 


LD-cavity,  ng  is  the  group  refraction  index  of  the  LD.).  Second,  fcxt  «  1  which  means 

only  one  round-trip  of  the  emitted  photons  is  significant  inside  the  EC.  Then  we  arrive  at 
the  following  two  rate  equations : 


m)  _ 

dr  ] 


a 


^s-^  +  ~gAN-Ns) 


£(0  +  ~(g--— )£(0  +  *£(f-r)exp(-/Qr) 


(1) 


dN  _  N  ' 

—  =  J - gS 

dt  T,P 


(2) 


where  g  —  g(N)(  1  fcpP)  is  the  nonlinear  gain  where  the  differential  gain  is  given  by 


g'(jV)  =  g N  { N  —  Nx ) ,  P  is  the  output  power,  Qv  is  the  resonance  optical  frequency  of 
the  LD  without  EC,  cc  is  the  linewidth  enhancement  factor,  g  v  is  the  differential  gain, 
tp  is  the  photon  lifetime  in  the  LD  cavity,  xsp  is  the  carrier  lifetime  due  to  spontaneous 
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and  nonradiative  recombination,  S  =  |£(r)|2  is  the  photon  number  in  the  LD  cavity, 
(output  power  P  S),  J  is  the  carrier  injection  rate.  The  feedback  parameter  is  given  by 


*0 


where  parameter  C,  which  is  a  function  of  the  facet  and  external  mirror  reflectivity  R2 
and  R2 


C, 


(1  -*2) 


is  a  measure  of  the  coupling  strength  between  the  two  cavities.  The  rate  equations  for 
the  EC  LD  are  given  in  equations  (1)  and  (2).  The  absolute  squared  value  of  the  electric 
field  amplitude  corresponds  to  the  number  of  photons  in  the  laser  cavity.  Equations  (1) 
and  (2)  can  be  cast  into  rate  equations  for  the  photon  number  S,  the  phase  (f>  ,  and  the 
carrier  number  N: 


dm 

dt 


)S(t)  +  2Ky[m)^fsm  r)cos (j>(t)  -  <j>(t  -  r)) 


(3) 


dm 

dt 


~  Ns)  —  rc 


4m 


sin  (nsT+</>(t)-0(t-r)) 


(4) 


dNjt) 

dt 


-gS(t) 


(5) 


3.2  Numerical  Simulations 

In  numerical  simulations,  the  system  of  rate  equations  system  (3) -(5)  is  solved  by 
using  a  fifth-order  Runge-Kutta  method  [3], [4].  We  specify  the  initial  conditions 
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S(0)  =  10,  ^(0)=0  and  N(0)=0.  In  the  time  interval  0<t<r,  the  data  of  photon 
number  S  and  phase  <f>  are  stored  in  memory  and  they  are  reused  after  the  round-trip 
time  r  of  the  EC  in  order  to  account  for  the  optical  feedback. 


Time  (p  s ) 


Fig.  6  The  strongest  oscillation  for  the  Lext=5.25mm  and  the  fal  group.  The 
feedback  parameter  is  0.206.  The  oscillation  frequency  is  equal  to  21.3  GHz. 


3.3  An  Example  of  High-Frequency  Microwave  Oscillation 

v 

For  example,  we  set  the  output  power  P=5  mW/per  facet.  The  calculations  used  the 
laser  parameters  given  in  Table  I.  Here,  we  also  choose  the  EC  length  Lex,  =5.25  mm. 
We  found  that  the  feedback  parameters  for  the  strongest  modulation  index  occurs  when 
/«,  =  0.206.  Further  increasing  fexl  destroys  the  twin  mode  locking.  It  leads  the  laser 
operation  on  a  new  single-mode  state  with  higher  output  power  (7.06  mW).  The 
strongest  oscillation  frequency  we  get  in  the  group  is  21.3  GHz  as  shown  in  Fig.  2. 
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IABLE  I 

List  of  Parameters  for  Fabry-Perot  laser 


Parameter 

Symbol 

Value 

Units 

scattering  loss  in  the  active  region 

3.6  x  103 

m-1 

carrier  number  at  transparency 

K 

2.0  x  10s 

output  power 

P 

2.906  x  10"8  x  S 

mW  !  facet 

index  of  active  region 

nx 

4.5 

differential  gain 

Sn 

2.2  xlO3 

s-x 

linewidth  enhancement  factor 

a 

4 

laser  wavelength 

X 

.3 

urn 

spontaneous  emission  lifetime 

r. 

2.0  xlO-9 

s 

photon  lifetime 

TP 

.29072  x  10~12 

s 

gain  saturation  coefficient 

K 

4.8 

w 1 

4.  Investigation  of  Deterministic  Properties:  Bifurcation 
Analysis 


Since  the  problem  is  nonlinear  in  nature,  there  does  not  exist  a  general  analytical 
methodology  that  can  lead  to  a  global  conclusion  of  the  operation  characteristics  of  the 
ECLD.  Our  approach  is  to  establish  a  numerical  technique  that  is  capable  of  predicting 
the  nature  of  the  ECLD  output  given  a  complete  set  of  parameters  for  the  ECLD 
without  numerical  integration  described  in  Section  2.  We  used  a  public-domain  software 
called  AUTO  to  investigate  the  dynamics  of  laser  diodes  with  short  external  cavity. 
AUTO  was  developed  in  the  Applied  Mathematics  Department  at  the  California 
Institute  of  Technology  as  one  of  the  projects  of  the  Differential  Equations  Group  in  the 
Center  for  Research  on  Parallel  Computation.  In  what  follows,  we  describe  the  interface 
FORTRAN  code  for  AUTO  to  simulate  the  laser  diodes  with  external  cavity  and  some 
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representative  results.  We  stress  that  the  results  shown  here  only  apply  to  the  particular 
parameters  chosen  as  indicated.  To  find  the  dynamic  behavior  of  ECLDs  with  different 
parameters,  one  needs  to  run  the  AUTO  program  with  the  interface  FORTRAN  code 
with  the  modified  ECLD  parameters. 

4.1  Bifurcation  Diagrams  of  ECLD  Steady  States 


To  construct  the  bifurcation  diagrams  of  ECLD  steady  states,  the  first  step  is  to  find  all 
ECLD  steady  states  for  a  given  set  of  the  ECLD  parameters  and  then  using  AUTO  to 
finding  the  continuations  of  the  steady  states. 

To  find  all  ECLD  steady  states  for  a  given  set  of  the  ECLD  parameters,  a 
program  such  as  SEARCH_SS  can  be  used.  SEARCH_SS  finds  all  the  solutions  to 
the  nonlinear  algebraic  equations: 

( g  — —)S (0  +  2 cos(Q, s t + -  <f>(t  -  r))  =  0  (6) 

Tp 

^-  —  -^5(0  =  0  (7) 

Tsp 

f 

arrived  by  setting  the  right-hand  side  of  Eqs.  (3)  and  (5)  to  be  equal  to  zero.  The 
program  “SEARCHES S”  is  given  in  Appendix  1.. 

Once  all  the  steady  states  of  the  ECLD  for  the  given  set  of  parameters  are  found, 
the  bifurcation  diagram  can  be  constructed  by  continuing  the  steady  states  as  one  of  the 
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parameters  in  the  parameter  set  is  varied  in  both  directions  (increasing  or  decreasing). 
The  interface  FORTRAN  code  for  this  continuation  is  contained  in  Appendix  2.  It  is 
found  that  Hopf  bifurcation  delineates  stable  and  unstable  steady  states. 

4.2  Stability  Boundaries  of  ECLD  Modes 

To  construct  the  Stability-Boundary  diagram  of  ECLD  steady  states,  the  first  step  is  to 
find  the  Hopf  bifurcation  points  of  all  ECLD  steady  states  for  a  given  set  of  the  ECLD 
parameters  and  then  using  AUTO  to  finding  the  continuations  of  the  Hopf  bifurcation 
points. 

To  find  Hopf  bifurcation  points  of  all  ECLD  steady  states  for  a  given  set  of  the 
ECLD  parameters,  a  program  such  as  “SEARCH_ST”  can  be  used.  SEARCH_ST  finds 
all  the  solutions  to  the  nonlinear  algebraic  equations  (6),  (7),  and  (8) 

V 

D(s)  =  - s 3  +  s2[y  +  2/ce(r )  cos^]-.sj<w  R 2  +tc2e 2  (r)  +  2ye/ce(r)cos^j 
2 e{z){oos(j)- asm <j))  +  y eKe2{z)  + ay j pe{r) sin (f\^ 

where  e(r)  =  l  - e'aT  and  <f>  =  Q.r.  Here  ye  =  1/  rs+gNS  is  the  carrier  perturbation 
decay  rate,  y  =  ye  +  yp  is  the  decay  rate  of  relaxation  oscillations,  with  y p  -  gKpP 
being  the  contribution  of  the  nonlinear  gain.  The  angular  frequency  of  relaxation 
oscillations  a>R  is  given  by  coR2  =  ggNS  +  yeyp.  Equation  (8)  is  arrived  by  setting  the 

determinant  of  the  linearized  ECLD  system  to  be  equal  to  zero.  The  program 
“SEARCH_SS”  is  given  in  Appendix  3. 
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SQRT  [rexr] 


Fig.  7  The  stability  boundaries  of  4  ECLD  modes  in  the  -J  fext  —  t d  plane.  The  injection 

current  is  at  twice  the  threshold. 

Once  Hopf  bifurcation  points  of  all  the  steady  states  of  the  ECLD  for  the  given 
set  of  parameters  are  found,  the  bifurcation  diagram  can  be  constructed  by  continuing 
the  steady  states  when  one  of  the  parameters  in  the  parameter  set  is  varied  in  both 
directions  (increasing  or  decreasing).  The  interface  FORTRAN  code  for  this 
continuation  is  contained  in  Appendix  4. 
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4.3  Regions  of  Existence  of  ECLD  Modes 
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Fig.  8  Limiting  points  (LP)  of  existence  of  modes  superimposed  on  the  stability 
boundaries.  When  LP  of  a  higher-order  mode  is  above  the  stability  boundary,  robust 
high-frequency  oscillations  result.  Parameters  are  the  same  as  in  Fig.  3,  except  that 

injection  is  three  time  the  threshold. 


It  is  found  by  numerical  simulation  that  even  in  regions  above  the  stability 
boundary  of  a  particular,  the  high-frequency  oscillations  are  not  always  present.  We 
believe  that  this  due  to  the  existence  of  other  modes,  which  might  be  stable  in  the 
unstable  region  of  the  mode  that  we  are  interested  in.  To  do  so,  we  calculated  the 
regions  of  existence  of  all  modes.  High-frequency  oscillations  are  produced  only  if  the 
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boundary  of  existence  [limiting  point  (LP)  of  existence]  of  higher-order  modes  is  above 
the  stability  boundary.  An  example  where  robust  high-frequency  oscillation  are  possible 
is  shown  below. 


4.4  Investigation  of  Deterministic  Properties:  Experiments 


-20  ; 
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*45 


■50 


Frequency 


Fig.  9  RF  spectrum  of  a  9.2  GFIz  subcarrier  generated  in  ECLD. 

Using  a  Fabry-Perot  laser  diode  and  an  optical  flat  as  the  external  mirror,  The  high- 
frequency  intensity  oscillations  have  been  generated.  The  highest  frequency  achieved  is 
9.5  GHz,  as  shown  in  Fig.  9.  Modulation  of  these  oscillations  has  also  been 
demonstrated.  Up  to  750  MHz  modulation  has  demonstrated.  The  modulated  signal, 
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shown  in  Fig.  10  reveals  excellent  frequency  modulation  characteristics  with  minimal 
residual  amplitude  modulation,  as  the  two  sidebands  are  almost  equal  in  strength. 


'Fig.  10  RF  spectrum  of  a  8.4  GHz  subcarrier  modulated  by  a  750  MHz  signal. 
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5.  Conclusions 


In  conclusion,  we  have  determine  the  underlying  mechanism  of  the  high-frequency 
intensity  oscillations  in  laser  diodes  with  short  external  cavity.  From  a  nonlinear  dynamic 
point  of  view,  Hopf  bifurcation  is  responsible  for  the  occurrence  of  the  high-frequency 
intensity  oscillations.  However,  Hopf  bifurcation  is  only  the  necessary  condition.  The 
necessary  and  sufficient  condition  are  that  Hopf  bifurcation  occurs  and  the  limiting  point 
of  existence  of  a  higher-order  mode  is  above  the  stability  boundary  of  a  lower-order 
mode.  The  noise  characteristics  of  the  high-frequency  oscillations  can  be  measured 
using  the  HP  3048A  phase  and  amplitude  noise  measurement  system. 
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Appendix  1 


PROGRAM  SEARCH  SS 


ECLD) . 


Search  for  steady  states  of  the 


external  cavity  laser  diode 


Program  based  on  algorithm  554  (routine  BRENTM)  of  TOMS. 


IMPLICIT  REAL* 8  (A-H,0-Z) 

DIMENSION  X(10) , FVEC(IO) ,WA(130) 

EXTERNAL  FUNC 

COMMON  /ECLD/  GN, TP, EP, R, TO, FEXT, TD, ALP, WO , TS, RJTH,  RJ,  WOT 
COMMON  /SEARCH/  ISEED, NSF, XF ( 100, 10 ), PAR ( 10) 

Get  ECLD  parameters  and  store  in  common  block  ECLD 

CALL  GETPAR 


The  C  parameter.  For  C  >>  1,  the  number  of  steady  states  is 
approximately  C/PI+1. 


WRITE  (*,100)  (1 . 0-R) *DSQRT (FEXT/R) *TD/T0 

OPEN  ( 1, FILE= ' search_ss . out 1 , STATUS= ' UNKNOWN ' ) 
WRITE  (1,*)  '  • 

CLOSE  (1) 


N  =  2 
LWA  =  130 


(number  of  equations 

(size  of  array  WA  (must  be  at  least  N*(N+3)) 


WRITE  (*,*)  'Enter  TOL, NTRY  (1 . OE-15, 1000  typical):' 
READ  (*,*)  TOL, NTRY 


ISEED  =  1  (seed  for  random  number  generator  in  START 

v  NSF  =  0  (number  of  solutions  found 

DO  10  I  =  1, NTRY 


Generate  trial  starting  values. 
CALL  START (N,X) 


Attempt  to  locate  a  root. 


CALL  BRENT1 (FUNC,N,X, FVEC, TOL, INFO, WA, LWA) 

If  current  search  is  successfull,  update  library  of  solutions. 
IF  ( INFO . LE . 3 . AND . INFO . NE . 0 )  CALL  UPDATE (N, X, INFO) 


10  CONTINUE 


Store  and  write  roots  in  ascending  order. 

CALL  PIKSRT (NSF,XF(1,1) , XF (1,2) ) 

OPEN  ( 1 , FILE= ' search_ss . out ' , STATUS= ' UNKNOWN ' ) 
DO  20  I  =  1 , NSF 
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Determine  if  the  current  solution  is  a  mode  or  antimode 
ET  =  ( 1 . 0-R) *DSQRT ( FEXT/R) *TD/T0 

DW  =  1 . 0+ET*DSQRT ( 1 . 0+ALP**2 ) *DCOS (XF (1,1) +DATAN (ALP) ) 

WRITE  (1,200)  XF (1,1), XF (1,2) , DW/DABS ( DW) 

20  CONTINUE 
CLOSE  (1) 

100  FORMAT  (’C  =',E11.4) 

200  FORMAT  (10E16.8) 

STOP 

END 

SUBROUTINE  GETPAR 


Get  ECLD  parameters  and  store  in  common  block  ECLD. 

IMPLICIT  REAL* 8  (A-H,0-Z) 

COMMON  /ECLD/  GN, TP, EP, R, TO, FEXT, TD, ALP, W0, TS, RJTH,  RJ,  WOT 
C 

OPEN  ( 1 ,  FILE=  f  ecld. in 1 , STATUS=  *  OLD ' ) 

READ  (1,*)  GN  ! gain  coefficient  (8.1E-19  cm**3/ps) 

READ  (1,*)  TP  ! photon  lifetime  (2  ps ) 

READ  (1,*)  EP  !nonlinear  gain  coefficient  (9.0E-18  cm**3) 

READ  (1,*)  R  !laser  facet  reflectivity  (0.31) 

READ  (1,*)  TO  ! solitary  laser  cavity  round  trip  time  (8  ps) 

READ  (1,*)  FEXT  ! feedback  fraction  (1.0E-04) 

READ  (1,*)  TD  Idelay  in  external  cavity  (20ps  for  3mm  EC) 

READ  (1,*)  ALP  ilinewidth  enhancement  factor  (6) 

READ  (1,*)  W0  ! f req  of  solitary  laser  (1.45E+03  rad/ps  for 

WL=1300  nm) 

READ  (1,*)  TS  ! carrier  lifetime  (2ns) 

READ  (1,*)  RJTH  ! current  density  at  threshold  (8.55E+14  cm**3/ps) 
READ  (1,*)  RJ  ! current  density  relative  to  threshold  (1.2) 

CLOSE  (1) 

C  v 

WOT  =  W0*TD-DMOD (W0*TD, 8 . 0 DO* DAT AN ( 1 . 0D0) ) 

FEXT  =  FEXT* *2 
C 

RETURN 

END 

C 

SUBROUTINE  FUNC (N, X, F, M) 

C  - 

c 

C  Solve  for  steady  state  frequency  of  the  external  cavity  laser 

C  diode  (ECLD).  See  Ritter  and  Haug,  JOSA  B,  Vol.  10,  p.  130,  1993. 

C 

C  X ( 1 )  =  W*TD 

C  X ( 2 )  -  ES/E0 

C 

C  M  is  the  component  of  F  to  be  returned,  N  is  the  number  of 

equations 

C  and  X  is  an  N-vector  containing  values  of  the  independent 

variables. 

C 

IMPLICIT  REAL* 8  (A-H,0-Z) 
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DIMENSION  X(N),F(N) 

COMMON  /ECLD/  GN,  TP,  EP, R, TO, FEXT, TD, ALP, WO,  TS,  RJTH,  RJ,  WOT 
C 

ET  =  (1. 0-R) *DSQRT (FEXT/R) *TD/T0 

EO  =  TP* (RJ-1.0) *RJTH 

GS  =  1 . O/TP-2 . 0*ET/TD*DCOS (X ( 1) ) 

DN  =  (GS/ ( 1 . 0-EP*X (2 ) *E0 ) -1 . 0/TP) /GN 
C 

F ( 1 )  =  X ( 1) - (W0*TD-W0T) -ALP*GN*DN*TD/2 . 0+ET*DSIN (X (1) ) 
F(2)  =  GS*X (2 ) *E0- (RJ-1 . 0 ) *RJTH+DN/TS 
C 

RETURN 

END 

C 

SUBROUTINE  START (N,X) 


Generate 
quations . 


random  starting  values 


for  X. 


N  is  the  number  of 


IMPLICIT  REAL* 8  (A-H,0-Z) 

COMMON  /SEARCH/  ISEED, NSF,XF(100, 10) , PAR(IO) 

DIMENSION  X (N) 

Generate  a  random  number  between  X1L  and  X1H  within  which  X(l)  is 
expected  to  lie.  RAND  is  in  the  range  (0,1). 

ISEED  =  MOD ( (2**10+3) *ISEED, 2**20) 

RAND  =  DBLE (ISEED) /2**20 

X1L  =  -4 . 0D1*DATAN ( 1 . 0D0 ) 

X1H  =  4 . 0D1*DATAN ( 1 . 0D0 ) 

X(l)  =  (X1H-X1L) *RAND+X1L 

Generate  another  random  number. 

ISEED  =  MOD ( ( 2**10+3 ) * ISEED,  2**20 ) 

RAND  =  DBLE (ISEED) /2**20 

X1L  =  0.0D0 

*  X1H  =  2.0D0 

X(l)  =  (X1H-X1L) *RAND+X1L 

RETURN 

END 

SUBROUTINE  UPDATE (N,X, INFO) 


Update  the  library  of  solutions. 

IMPLICIT  REAL *8  (A-H,0-Z) 

DIMENSION  X (N) , XOLD (10) 

COMMON  /SEARCH/  ISEED, NSF, XF ( 100, 10) , PAR ( 10) 

Check  to  see  if  the  current  data  point  already  found  earlier  to 
within  a  relative  tolerance  of  1.0E-07.  If  so,  return. 

OPEN  ( 1 , FILE= ' search_ss . out ' , STATUS= ' OLD '  ) 

REWIND  1 

10  READ  (1, 100, END=30)  (XOLD(I),  1=1, N) 

DO  20  I  =  1 , N 

DIFF  =  DABS (XOLD(I) -X (I) ) 

IF  ( DIFF . GT .  ( 1 . 0D-07  *DABS (X ( I ) ) ) )  THEN 
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IF  (DABS (X (I) ) .GT.1.0D-10)  GO  TO  10 
ENDIF 
20  CONTINUE 
CLOSE  (1) 

RETURN 

C 

30  NSF  =  NSF+1 
C 

DO  40  I  =  1,N 

XF  (NSF,  I)  =  X  (I ) 

40  CONTINUE 
C 

OPEN  (1, FILE= ' search_ss . out  * , STATUS^ 1  OLD ’  ) 
WRITE  (1,100)  (X ( I ) ,  1=1, N) 

CLOSE  (1) 

C 

100  FORMAT ( 125E1 6 . 8 ) 

C 

RETURN 

END 


C 

C 

c 

c 

c 

of 

C 

C 


c 


c 


SUBROUTINE  PIKSRT  (N,ARR,BRR) 


Sorts  an  array  ARR  of  length  N  into  ascending  numerical  order,  by 
straight  insertion,  while  making  the  corresponding  rearrangement 

the  array  BRR. 

IMPLICIT  REAL *8  (A-H,0-Z) 

DIMENSION  ARR (N) , BRR (N) 


DO  30  J  =  2 , N 
A  =  ARR(J) 

B  =  BRR(J) 

DO  10  I  ~  J-1,1,-1 

IF  (ARR (I) . LE. A)  GO  TO  20 
1  ARR  ( I  + 1 )  ==  ARR  ( I ) 

BRR ( 1+1 )  -  BRR (I) 


10 

CONTINUE 

1  =  0 

20 

ARR { I + 1 )  =  A 

BRR ( 1+1 )  =  B 

30  CONTINUE 


RETURN 

END 


26 


noon  oo  oooo  ooooo  o  o  oooo  mood 


Appendix  2 


SUBROUTINE  FUNC (N, X, ICP, PAR, IJAC, F,  DFDU,  DFDP) 


C 

C  Continuation  of  the  steady  state  solution  of  the  external  cavitv 
C  laser  diode  (ECLD) . 

C 

C  X ( 1 ) =WT ,  where  WT  is  the  steady  state  phase  shift  in  external 

cavity 

C  X(2)=ES/E0,  where  ES  and  EO  are  the  steady  state  field 

intensities  of 

the  ECLD  and  the  solitary  laser,  respectively. 

PAR ( 1 ) =RJ ,  current  density  relative  to  threshold 

MPLICIT  REAL* 8  (A, B, D-H, O-Z ) 

IMPLICIT  COMPLEX* 16  (C) 

DIMENSION  X(N) ,F(N) ,PAR(*) 

COMMON  /ECLD/  GN, TP, EP, R, TO, FEXT, TD, ALP, WO, TS, RJTH,  RJ,  WOT 
Update  ECLD  parameters. 

RJ  =  PAR ( 1 ) 

ET  =  (1.0-R) *DSQRT (FEXT/R) *TD/TO 

EO  =  TP* (RJ-1.0) *RJTH 

GS  =  1 . O/TP-2 . 0*ET/TD*DCOS (X ( 1 ) ) 

DN  =  (GS/ ( 1 . 0-EP*X (2 ) *E0) -1 . 0/TP) /GN 

F ( 1 )  =  X ( 1 ) - (W0*TD-W0T) -ALP*GN*DN*TD/2 . 0+ET*DSIN (X (1)  ) 

F(2)  =  GS*X(2) *E0- (RJ-1.0) *RJTH+DN/TS 

RETURN 

END 

SUBROUTINE  STPNT ( NDIM, U, PAR) 


Starting  conditions. 

IMPLICIT  REAL *8  (A-H,0-Z) 

DIMENSION  U (NDIM) , PAR (*) 

COMMON  /ECLD/  GN, TP, EP, R, TO, FEXT, TD, ALP, WO, TS, RJTK, RJ,  WOT 
Get  ECLD  parameters  and  store  in  common  block  ECLD 
CALL  GETPAR 

OPEN  ( 1 , FILE= ' search_ss . out ' , STATUS= ' OLD '  ) 

READ  (1,*)  (U(I),  1=1, NDIM) 

CLOSE  (1) 

PAR ( 1 )  =  RJ 

RETURN 

END 

SUBROUTINE  GETPAR 


Get  ECLD  parameters  and  store  in  common  block  ECLD. 
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IMPLICIT  REAL *8  (A-H,0-Z) 

COMMON  /ECLD/  GN, TP, EP, R, TO, FEXT, TD, ALP, WO , TS,  RJTH,  RJ,  WOT 


gain  coefficient  (8.1E-19  cm**3/ps) 
photon  lifetime  (2  ps) 

nonlinear  gain  coefficient  (9.0E-18  cm**3) 

laser  facet  reflectivity  (0.31) 

solitary  laser  cavity  round  trip  time  (8  ps) 

feedback  fraction  (1.0E-04) 

delay  in  external  cavity  (20ps  for  3mm  EC) 

linewidth  enhancement  factor  (6) 

!  freq  of  solitary  laser  (1.45E+15  rad/s 


OPEN 

(i, 

FILE=  *  e 

READ 

(i, 

*) 

GN 

READ 

(i, 

*) 

TP 

READ 

(i, 

*> 

EP 

READ 

(i, 

*) 

R 

READ 

(i. 

*) 

TO 

READ 

(i, 

*) 

FEXT 

READ 

(i, 

*) 

TD 

READ 

a, 

*) 

ALP 

READ 

(i 

,*) 

WO 

WL=1300  nm) 

READ 

(i. 

*> 

TS 

READ 

(i. 

*) 

RJTH 

READ 

(i. 

*) 

RJ 

for 


CLOSE  (1) 


! carrier  lifetime  (2ns) 

! current  density  at  threshold  (8.55E+14  cm**3/ps) 
! current  density  relative  to  threshold  (1.2) 


WOT  =  WO*TD-DMOD (W0*TD, 8 . 0D0*DATAN (1 . 0D0)  ) 
FEXT  =  FEXT* *2 

RETURN 

END 


SUBROUTINE  BOND 


Dummy  routine 

RETURN 

END 

SUBROUTINE  ICND 


Dummy  routine 

RETURN 

END 

SUBROUTINE  FOPT 


Dummy  routine 

RETURN 

END 
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’  DW  -  1 . 0+ET*DSQRT ( 1 . O+ALP*  *2 ) *DCOS (XF (1,3) +DATAN (ALP) ) 

WRITE  (1,100)  XF {1,1) , XF (1,2) , XF (1,3) , XF (1,4),  DW/DABS ( DW) 
20  CONTINUE 
CLOSE  (1) 

C 

100  FORMAT  (11E16.8) 

C 

STOP 

END 

C 

SUBROUTINE  GETPAR 
C  - - 

c 

c  Get  ECLD  parameters  and  store  in  common  block  ECLD. 

C 

IMPLICIT  REAL* 8  (A-H,0-Z) 

COMMON  /ECLD/  GN, TP, EP, R, TO, FEXT, TD, ALP, W0,  TS,  RJTH,  RJ,  WOT 


OPEN  ( 1 , FILE= 1 ecld .in’, STATU S= ' OLD ' ) 

READ  (1,*)  GN  ! gain  coefficient  (8.1E-19  cm**3/ps) 

Iphoton  lifetime  (2  ps) 

! nonlinear  gain  coefficient  (9.0E-18  cm**3) 

! laser  facet  reflectivity  (0.31) 

! solitary  laser  cavity  round  trip  time  (8  ps) 

READ  (1,*)  FEXT  ! feedback  fraction  (1.0E-Q4) 

READ  (1,*)  TD  ! delay  in  external  cavity  (20ps  for  3mm  EC) 

READ  (1,*)  ALP  llinewidth  enhancement  factor  (6) 

Ifreq  of  solitary  laser  (1.45E+15  rad /s  for 


READ  (1,*)  TP 
READ  (1,*)  EP 
READ  ( 1 , * )  R 
READ  (1,*)  TO 


READ  ( 1 ,  *  )  W0 
WL=1300  nrn) 

READ  ( 1 , * )  TS 


! carrier  lifetime  (2ns) 

READ  (1,*)  RJTH  ! current  density  at  threshold  (8.55E+14  cm**3/ps) 
READ  (1,*)  RJ  ! current  density  relative  to  threshold  (1.2) 

CLOSE  (1) 

WOT  =  W0*TD-DMOD (W0*TD, 8 . 0 DO* DAT AN { 1 . 0D0 ) ) 

FEXT  =  FEXT* *2 

RETURN 

END 


C 
C 
c 
c 

(ECLD) 

C 

c 

c 

oscil . 
C 

cavity 

C 


SUBROUTINE  FUNC (N, X, F, M) 

Solve  for  steady  state  solutions  aand  roots  of  determinant  of 

linearized  rate  equations  of  the  external  cavity  laser  diode 

X ( 1 ) =SQRT ( FEXT ) ,  where  FEXT  is  the  feedback  fraction 
X ( 2 ) -W/2/PI ,  where  W  is  the  (circular)  frequency  of  relaxation 

X(3)— WT,  where  WT  is  the  steady  state  phase  shift  in  external 

X(4)—ES/E0,  where  ES  and-  E0  are  the  steady  state  field 

the  ECLD  and  the  solitary  laser,  respectively. 


intensities  of 
C 
C 

C  M  is  the  component  of  F  to  be  returned,  N  is  the  number  of 

equations 

C  and  X  is  an  N-vector  containing  values  of  the  independent 

variables . 
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Appendix  3 


PROGRAM  SEARCH  ST 

C  - 

c 

C  Search  for  stability  boundary  of  the  external  cavity  laser  diode 

( ECLD)  . 

C 

C  Program  based  on  algorithm  554  (routine  BRENTM)  of  TOMS. 

IMPLICIT  REAL *8  (A-H,0-Z) 

DIMENSION  X (10) , FVEC{10) ,WA(130) 

EXTERNAL  FUNC 

COMMON  /ECLD/  GN, TP, EP, R, TO, FEXT, TD, ALP, WO,  TS,  RJTH,  RJ,  WOT 
COMMON  /ECLDSS/  WS,  ES,  DN, GAMMAN, GAMMAI , GAMMA, OMEGSQ, ALPT 
COMMON  /SEARCH/  ISEED, NSF, XF ( 500 , 10) , PAR (10) 


Get  ECLD  parameters  and  store  in  common  block  ECLD 


C 


C 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


CALL  GET PAR 


OPEN  ( 1 ,  FILE- 1 search_s t . out 1 , STATUS^ 1  UNKNOWN ' ) 
WRITE  (1,+)  »  1 

CLOSE  (1) 


N  =  4 
LWA  -  130 


!  number  of  equations 

! size  of  array  WA  (must  be  at  least  N*(N+3)) 


WRITE  (*,*)  'Enter  TOL, NTRY  ( 1 . OE-15,  1000  typical ): 1 
READ  (*,*)  TOL, NTRY 


ISEED  -  1 
NSF  =  0 

DO  10  I  =  1 , NTRY 


!seed  for  random  number  generator  in  START 
!  number  of  solutions  found  so  far 


Generate  trial  starting  values. 

CALL  START (N,X) 

Attempt  to  locate  a  root. 

CALL  BRENT1 ( FUNC, N, X, FVEC, TOL, INFO, WA, LWA) 

If  current  search  is  successfull,  update  library  of  solutions. 

IF  ( INFO . LE . 3 . AND . INFO . NE . 0 )  THEN 

IF  (X  (1 )  . GT. 0. 0D0.AND.X (2)  . GT . 1 . 0D-04 . AND . X ( 4 )  .GT.O.ODO 
.AND.X(l) .LT.0.5D0)  CALL  UPDATE (N, X, INFO) 

ENDIF 


10  CONTINUE 


Sort  and  write  roots  in  ascending  order. 

CALL  PIKSRT (NSF,XF(1, 1 ) , XF ( 1 , 2 ) , XF ( 1 , 3 ) , XF ( 1 , 4 ) ) 

OPEN  (1, FILE=T search_st . out ' , STATUS^ ' UNKNOWN ’ ) 

DO  20  I  ~  1 , NSF 

Determine  if  the  current  solution  is  a  mode  or  antimode 


ET  -  (1.0~R)*XF(I,1) /DSQRT (R) *TD/T0 
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c 


c 


IMPLICIT  REAL *8  (A, B, D-H, 0-Z) 

IMPLICIT  COMPLEX* 16  (C) 

DIMENSION  X (N) , F(N) 

COMMON  /ECLD/  GN, TP, EP, R, TO, FEXT, TD, ALP, WO, TS, RJTH,  RJ,  WOT 
COMMON  /ECLDSS/  WS, ES , DN, GAMMAN, GAMMAI , GAMMA, OMEGSQ, ALPT 

FEXT  =  X (1) **2 


Determine  steady  state  parameters  and  store  in  common  block  ECLDSS 
WS  =  X ( 3 ) 

EO  =  TP* (RJ-1.0) *RJTH 
ES  =  X ( 4 ) *E0 

ET  =  ( 1 . 0-R) *X ( 1 ) /DSQRT (R) *TD/T0 
GS  =  1 . O/TP-2 . 0*ET/TD*DCOS (WS) 

DN  =  (GS/ (1 . 0-EP*ES) -1.0/TP) /GN 
CALL  STDST 


Complex  frequency  and  coupling  strength 

CZ  =  DCMPLX { 0 . ODO, 8 . 0D0*DATAN ( 1 . 0D0) *X (2 ) ) 
CE  =  ET/TD* (1 . 0-CDEXP (-CZ*TD) ) 


The  complex  determinant 
CDET  =  CZ**3+ 

*  2.0* ( GAMMA+CE * DCOS (WS) ) *CZ**2+ 

*  (OMEGSQ+2 . 0* ( GAMMA+GAMMAN ) *CE*DCOS (WS) +CE**2) *CZ+ 

*  2 . 0*GAMMAN*CE**2+ 

*  OMEGSQ* DSQRT (1 . 0+ALPT**2) *CE*DCOS (WS+DATAN (ALPT) ) 

F ( 1 )  =  DREAL (CDET) 

F (2)  =  DIMAG (CDET) 

F(3)  =  WS- (W0*TD-W0T) -ALP*GN*DN*TD/2 . 0+ET*DSIN (WS) 

F ( 4 )  =  GS*ES- (RJ-1.0) *RJTH+DN/TS 

RETURN 
v  END 

SUBROUTINE  STDST 


Determine  steady  state  parameters  of  ECLD.  The  steady  state  data 
stored  in  common  block  ECLDSS. 

IMPLICIT  REAL *8  (A-H,0-Z) 

COMMON  /ECLD/  GN, TP, EP, R, TO, FEXT, TD, ALP, WO , TS, RJTH,  RJ,  WOT 
COMMON  /ECLDSS/  WS, ES, DN, GAMMAN, GAMMAI, GAMMA,  OMEGSQ,  ALPT 

Damping  GAMMAN  caused  by  spontaneous  and  induced  emission 

GNP  =  GN* ( 1 . 0-EP*ES) 

GAMMAN  =  (1.0+GNP*ES*TS) /2.0/TS 

Damping  GAMMAI  caused  by  nonlinear  gain  saturation 

EPP  =  EP/ ( 1 . 0-EP*ES ) 

GS  =  (1. 0-EP*ES) * (GN*DN+1 . 0/TP) 

GAMMAI  =  EPP*ES*GS/2 . 0 
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Decay  rate  GAMMA  of  relaxation  oscillations  of  solitary  laser 
GAMMA  =  G AMMAN + GAMMAI 

Square  of  relaxation  oscillation  frequency  0MEG**2 

OMEGSQ  =  GNP*GS*ES* (1 . O+EPP/GNP/TS) 

(Modified)  linewidth  enhancement  factor. 

ALPT  =  ALP/ (1 . 0-EP*ES) * (1 . 0-4 . 0*GAMMAN*GAMMAI/OMEGSQ) 

RETURN 

END 

SUBROUTINE  START (N,X) 


Generate  random  starting  values  for  X.  N  is  the  number  of 
tions. 

IMPLICIT  REAL *8  (A-H,0-Z) 

COMMON  /SEARCH/  ISEED, NSF, XF ( 500, 10) , PAR ( 10) 

DIMENSION  X (N) 

Generate  a  random  number  between  X1L  and  X1H  within  which  X(l)  is 
expected  to  lie.  RAND  is  in  the  range  (0,1). 

ISEED  =  MOD{ (2 **10+3) * I SEED,  2**20) 

RAND  =  DBLE (ISEED) /2**20 

X1L  =  0.0D+00 

X1H  =  1.0D+00 

X(l)  =  (X1H-X1L) *RAND+X1L 

Another  random  number  for  X(2). 

ISEED  =  MOD( (2** 10+3) * ISEED,  2**20) 

•  RAND  =  DBLE (ISEED) /2**20 
X2L  =  1.0D-04 

X2H  =  1.0D-01 

X (2)  =  (X2H-X2L) *RAND+X2L 

Another  random  number  for  X(3). 

ISEED  =  MOD ( { 2  *  * 10+3 ) *ISEED,  2**20) 

RAND  =  DBLE (ISEED) /2**20 
X3L  =  -4 . 0D1*DATAN ( 1 . 0D0 ) 

X3H  =  4 . 0D1*DATAN (1 . 0D0) 

X (3)  =  (X3H-X3L) *RAND+X3L 

Another  random  number  for  X(4). 

ISEED  =  MOD ( (2** 10+3) * ISEED,  2**20) 

RAND  =  DBLE (ISEED) /2**20 
X4L  =  0.5D0 
X4H  =  2.5D0 

X ( 4 )  =  (X4H-X4L) *RAND+X4L 

RETURN 

END 


32 


oooonooo  oo  n  no  oooo  noon 


SUBROUTINE  UPDATE (N,  X, INFO) 


Update  the  library  of  solutions. 

IMPLICIT  REAL *8  (A-H,0-Z) 

DIMENSION  X(N) ,XOLD( 10) 

COMMON  /SEARCH/  ISEED, NSF, XF ( 500, 10) , PAR ( 10) 

Check  to  see  if  the  current  data  point  already  found  earlier  to 
within  a  relative  tolerance  of  1.0E-07.  If  so,  return. 

OPEN  (1, FILE=' search  st . out ' , STATUS= ' OLD ' ) 

REWIND  1 

10  READ  (1, 100,END=30)  (XOLD(I),  1=1, N) 

DO  20  I  =  1,N 

DIFF  =  DABS (XOLD ( I ) -X ( I ) ) 

IF  (DIFF. GT . (1 . 0D- 07* DABS (X ( I ) ) ) )  THEN 
IF  (DABS (X (I) ) . GT. 1 . 0D-10)  GO  TO  10 
ENDIF 
20  CONTINUE 
CLOSE  (1) 

RETURN 

30  NSF  =  NSF+1 

DO  40  I  =  1 , N 

XF  (NSF,  I )  =  X  (I) 

40  CONTINUE 

OPEN  (1, FILE= ' search  st. out' , STATUS='OLD' ) 

WRITE  (1,100)  (X (I) ,  1=1, N) 

CLOSE  (1) 

100  FORMAT (125E16. 8) 

RETURN 
*  END 

SUBROUTINE  PIKSRT  (N, ARR, BRR, CRR, DRR) 


Sorts  an  array  ARR  of  length  N  into  ascending  numerical  order,  by 
straight  insertion,  while  making  the  corresponding  rearrangement 

the  arrays  BRR,  CRR  and  DRR. 

IMPLICIT  REAL  *  8  (A-H,0-Z) 

DIMENSION  ARR(N) , BRR (N) ,CRR(N) , DRR (N) 

DO  30  J  =  2, N 
A  =  ARR(J) 

B  =  BRR(J) 

C  =  CRR(J) 

D  =  DRR ( J) 

DO  10  I  =  J-1,1,-1 

IF  (ARR ( I ) . LE . A)  GO  TO  20 
ARR ( I + 1 )  =  ARR (I) 

BRR (1+1 )  =  BRR ( I ) 

CRR ( 1+1 )  =  CRR ( I ) 
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DRR ( I + 1 ) 


10 

CONTINUE 
1  =  0 

20 

ARR { I + 1 ) 

=  A 

BRR { 1  +  1 ) 

=  B 

CRR ( 1+1 ) 

=  C 

DRR ( 1  +  1 ) 

=  D 

30 

CONTINUE 

RETURN 

END 

v 


DRR  ( I ) 
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Appendix  4 


SUBROUTINE  FUNC (N, X, ICP, PAR, IJAC, F, DFDU, DFDP) 

C  - - 

c 

C  Roots  of  determinant  of  linearized  rate  equations  of  the  external 
C  cavity  laser  diode  (ECLD) . 

C 

C  X (1) =SQRT (FEXT) ,  where  FEXT  is  the  feedback  fraction 

C  X(2)=W/2/PI,  where  W  is  the  (circular)  frequency  of  relaxation 

oscil. 

C  X ( 3 ) ~WT ,  where  WT  is  the  steady  state  phase  shift  in  external 

cavity 


C  X ( 4 ) =ES/E0 , 

intensities  of 

C  the  ECLD  and  the  solitary  laser,  respectively. 

C  PAR ( 1 ) -TD-TDO,  where  TD  is  round  trip  delay  in  external  cavity 

in  ps, 

C 
C 


where  ES  and  EO  are  the  steady  state  field 


and  TDO  is  starting  delay  in  ps. 


C 

C 

C 

C 


C 

c 

c 


c 

c 

c 


c 

c 

c 


c 


IMPLICIT  REAL *8  (A, B, D-H, 0-Z) 

IMPLICIT  COMPLEX* 16  (C) 

DIMENSION  X ( N ) ,  F  ( N ) , PAR ( * ) 

COMMON  /ECLD/  GN,  TP, EP, R, TO, FEXT, TD, ALP, WO, TS, RJTH, RJ,  WOT,  TDO 
COMMON  /ECLDSS/  WS ,  ES , DN, GAMMAN, GAMMAI , GAMMA, OMEGSQ, ALPT 

PI  =  4 . 0 DO* DAT AN ( 1 . 0D0 ) 

Update  ECLD  parameters. 

TD  =  PAR(1)+TD0  ! delay  in  external  cavity 

R J  =  PAR ( 1 ) 

FEXT  =  X(l)**2  [feedback  fraction 

Determine  steady  state  parameters  and  store  in  common  block  ECLDSS 
WS  =  X ( 3 ) 

EO  =  TP* (RJ-1 . 0) *RJTH 
'  ES  =  X ( 4 ) *E0 

ET  =  ( 1 . 0-R) *X ( 1 ) /DSQRT (R) *TD/T0 
GS  =  1 . O/TP-2 . 0*ET/TD*DCOS (WS) 

DN  =  (GS/ (1.0-EP*ES)-1.0/TP) /GN 
CALL  STDST 


Complex  frequency  and  coupling  strength 
CZ  =  DCMPLX  ( 0 . 0D0, 2 . 0D0*PI*X (2)  ) 

CE  =  ET/TD* (1 . 0-CDEXP (-CZ*TD) )  ^  : 

The  complex  determinant 
CDET  =  CZ**3+ 

2 . 0* ( GAMMA+CE*  DCOS (WS) ) *CZ**2  + 

(OMEGSQ+2 . 0* (GAMMA+GAMMAN) *CE*DCOS (WS) +CE**2) *CZ+ 
2. 0*GAMMAN*CE**2+ 

OMEGSQ*DSQRT (1 . 0+ALPT**2) *CE*DCOS (WS+DATAN (ALPT) ) 

F ( 1 )  =  DREAL (CDET) 

F ( 2 )  =  DIMAG (CDET) 

F ( 3 )  =  WS- (WQ*TD-W0T) -ALP*GN*DN*TD/2 . 0+ET*DSIN (WS) 

F ( 4 )  =  GS*ES- (RJ-1 . 0) *RJTH+DN/TS 


ooo  o  o  s;  o  ooono  o  o  - o  o  o  o  oono 


C 

RETURN 

END 

C 

SUBROUTINE  STPNT (NDIM, U, PAR) 


Starting  conditions . 

IMPLICIT  REAL* 8  (A-H,0-Z) 

DIMENSION  U (NDIM) , PAR (*) 

COMMON  /ECLD/  GN,  TP,  EP, R, TO, FEXT, TD, ALP, WO, TS, RJTH,  RJ,  WOT,  TDO 

Get  ECLD  parameters  and  store  in  common  block  ECLD 

CALL  GETPAR 

TDO  -  TD 

PAR ( 1 )  =  TD-TDO 

PAR ( 1 )  =  RJ 

OPEN  ( 1 , FILE= 1  searches t . out 1 , STATUS^ 1  OLD 1  ) 

READ  (1,*)  ( U ( I ) ,  1=1, NDIM) 

CLOSE  (1) 

RETURN 

END 

SUBROUTINE  GETPAR 


Get  ECLD  parameters  and  store  in  common  block  ECLD. 


IMPLICIT  REAL *8  (A-H,0-Z) 

COMMON  /ECLD/  GN, TP, EP, R, TO , FEXT, TD, ALP, WO , TS , RJTH, RJ,  WOT, TDO 


OPEN 

(i, 

FILE= 'ecld .in’, STATUS= ’ OLD ' ) 

READ 

(i, 

*) 

GN 

! gain  coefficient  (8. IE-19  cm**3/ps) 

'  READ 

(i, 

*) 

TP 

{photon  lifetime  (2  ps) 

READ 

(i, 

*) 

EP 

! nonlinear  gain  coefficient  (9.0E-18  cm 

**3) 

READ 

(i, 

*) 

R 

{laser  facet  reflectivity  (0.31) 

READ 

(i, 

*) 

TO 

! solitary  laser  cavity  round  trip  time 

(8  ps) 

READ 

(i, 

*) 

FEXT 

{feedback  fraction  (1.0E-04) 

READ 

(i, 

*) 

TD 

! delay  in  external  cavity  (20ps  for  3mm 

EC) 

READ 

(i. 

*) 

ALP 

llinewidth  enhancement  factor  (6) 

READ 

(i 

,*) 

WO 

!  f req  of  solitary  laser  (1.45E+15 

rad/s  for 

,=  1300  nm) 

READ 

(i, 

*) 

TS 

{carrier  lifetime  (2ns) 

READ 

a. 

*) 

RJTH 

{current  density  at  threshold  (8.55E+14 

cm**3/ps) 

READ 

(i, 

*) 

RJ 

{current  density  relative  to  threshold 

(1.2) 

CLOSE  (1) 

WOT  =  WO*TD-DMOD (W0*TD, 8 . ODO*DATAN ( 1 . 0D0 ) ) 
FEXT  =  FEXT* *2 

RETURN 

END 


SUBROUTINE  STDST 
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c 

are 

C 

C 


C 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


Determine  steady  state  parameters  of  ECLD.  The  steady  state  data 
stored  in  common  block  ECLDSS. 

IMPLICIT  REAL *8  (A-H,0-Z) 

COMMON  /ECLD/  GN, TP, EP, R, TO, FEXT, TD, ALP, WO,  TS, RJTH,  RJ,  WOT,  TDO 
COMMON  /ECLDSS/  WS,  ES, DN, GAMMAN, GAMMAI , GAMMA, OMEGSQ,  ALPT 

Damping  GAMMAN  caused  by  spontaneous  and  induced  emission 

GNP  =  GN* (1.0-EP*E$) 

GAMMAN  =  (1.0+GNP*ES*TS)/2.0/TS 

Damping  GAMMAI  caused  by  nonlinear  gain  saturation 

EPP  =  EP/ (1. 0-EP*ES) 

GS  =  (1.0-EP*ES)*(GN*DN+1.0/TP) 

GAMMAI  =  EPP*ES*GS/2 . 0 

Decay  rate  GAMMA  of  relaxation  oscillations  of  solitary  laser 
GAMMA  -  GAMMAN + GAMMAI 

Square  of  relaxation  oscillation  frequency  OMEG**2 

OMEGSQ  -  GNP*GS*ES* (1.0+EPP/GNP/TS) 

(Modified)  linewidth  enhancement  factor. 

ALPT  =  ALP/ ( 1 . 0-EP*ES ) * ( 1 . 0-4 . 0* GAMMAN* GAMMAI /OMEGSQ) 

RETURN 

END 

SUBROUTINE  BOND 

Dummy  routine 

RETURN 

END 

SUBROUTINE  ICND 

Dummy  routine 

RETURN 

END 

SUBROUTINE  FOPT 


Dummy  routine 

RETURN 

END 
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